#include <gdal/gdalwarper.h>
#include <gdal/cpl_conv.h>
#include <gdal/ogr_spatialref.h>

int main(int argc, char **argv)
{
	// 打开输入输出图像
	GDALAllRegister();

	GDALDriverH hDriver;
	GDALDataType eDT;
	GDALDatasetH hDstDS;
	GDALDatasetH hSrcDS;
	int nBand = 1;

	// 打开源文件
	hSrcDS = GDALOpen("N20E125.tif", GA_ReadOnly);
	CPLAssert(hSrcDS != NULL);

	// 创建输出图像的数据类型和输入图像第一个波段类型一致
	eDT = GDALGetRasterDataType(GDALGetRasterBand(hSrcDS, nBand));

	// 获取输出图像的驱动（GeoTIFF格式）
	hDriver = GDALGetDriverByName("GTiff");
	CPLAssert(hDriver != NULL);

	// 获取源图像的坐标系统
	const char *pszSrcWKT = GDALGetProjectionRef(hSrcDS);
	CPLAssert(pszSrcWKT != NULL && strlen(pszSrcWKT) > 0);

	// 设置输出图像的坐标系统为WGS84
	char *pszDstWKT = NULL;
	OGRSpatialReference oSRS;
	// oSRS.SetUTM(11, TRUE);
	// oSRS.SetWellKnownGeogCS("WGS84");
	oSRS.importFromEPSG(4326);
	oSRS.exportToWkt(&pszDstWKT);

	// 创建一个变换参数，从源图像的行列号坐标到结果图像的地理坐标（没有
	// 结果行列）的句柄，初始值设置为NULL。
	void *hTransformArg;
	hTransformArg = GDALCreateGenImgProjTransformer(hSrcDS, pszSrcWKT, NULL, pszDstWKT,
													FALSE, 0, 1);
	CPLAssert(hTransformArg != NULL);

	// 获取输出文件的近似地理范围和分辨率
	double adfDstGeoTransform[6];
	int nPixels = 0, nLines = 0;
	CPLErr eErr;
	eErr = GDALSuggestedWarpOutput(hSrcDS,
								   GDALGenImgProjTransform, hTransformArg,
								   adfDstGeoTransform, &nPixels, &nLines);
	CPLAssert(eErr == CE_None);

	GDALDestroyGenImgProjTransformer(hTransformArg);

	// 创建输出文件
	hDstDS = GDALCreate(hDriver, "out.tif", nPixels, nLines,
						GDALGetRasterCount(hSrcDS), eDT, NULL);

	CPLAssert(hDstDS != NULL);

	// 写入投影
	GDALSetProjection(hDstDS, pszDstWKT);
	GDALSetGeoTransform(hDstDS, adfDstGeoTransform);

	GDALRasterBandH hSrcRasterBand = GDALGetRasterBand(hSrcDS, nBand);
	GDALRasterBandH hDstRasterBand = GDALGetRasterBand(hDstDS, nBand);
	// 复制颜色表，如果有的话
	GDALColorTableH hCT;
	hCT = GDALGetRasterColorTable(hSrcRasterBand);
	if (hCT != NULL)
		GDALSetRasterColorTable(hDstRasterBand, hCT);
	GDALSetRasterNoDataValue(hDstRasterBand, GDALGetRasterNoDataValue(hSrcRasterBand, &nBand));

	// 建立变换选项
	GDALWarpOptions *psWarpOptions = GDALCreateWarpOptions();

	psWarpOptions->hSrcDS = hSrcDS;
	psWarpOptions->hDstDS = hDstDS;

	psWarpOptions->nBandCount = 1;
	psWarpOptions->panSrcBands =
		(int *)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
	psWarpOptions->panSrcBands[0] = nBand;
	psWarpOptions->panDstBands =
		(int *)CPLMalloc(sizeof(int) * psWarpOptions->nBandCount);
	psWarpOptions->panDstBands[0] = nBand;
	// psWarpOptions->papszWarpOptions = {"INIT_DEST=NO_DATA", ""};
	psWarpOptions->pfnProgress = GDALTermProgress;

	// 创建重投影变换函数
	psWarpOptions->pTransformerArg =
		GDALCreateGenImgProjTransformer(hSrcDS,
										GDALGetProjectionRef(hSrcDS),
										hDstDS,
										GDALGetProjectionRef(hDstDS),
										FALSE, 0.0, 1);
	psWarpOptions->pfnTransformer = GDALGenImgProjTransform;

	// 初始化并且执行变换操作
	GDALWarpOperation oOperation;

	oOperation.Initialize(psWarpOptions);
	oOperation.ChunkAndWarpImage(0, 0,
								 GDALGetRasterXSize(hDstDS),
								 GDALGetRasterYSize(hDstDS));

	GDALDestroyGenImgProjTransformer(psWarpOptions->pTransformerArg);
	GDALDestroyWarpOptions(psWarpOptions);

	GDALClose(hDstDS);
	GDALClose(hSrcDS);

	return 0;
}